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A novel view of plane wave expansion method in photonic crystals 

Young-Chung Hsu^] and Tzong-Jer Yanj^j] 
Department of Electrophysics, National Chiao-Tung University, Hsinchu, Taiwan, Republic of China 

(Dated: February 2, 2008) 

We propose a method derived from the simple plane wave expansion that can easily solve the 
interface problem between vacuum and a semi-infinite photonic crystal. The method is designed to 
find the complete set of all the eigenfunctions, propagating or evanescent, of the translation operators 
{Tr}, at a fixed frequency. With these eigenfunctions and their eigenvalues, the transmitted and 
reflected waves can be determined. Two kinds of applications are presented for 2D photonic crystals. 
The first is a selection rule for determine the normal direction of the vacuum-photonic crystal 
interface to achieve the highest attenuation effect at a gap frequency. The second is to calculate the 
transmittance and reflectance for a light incident from vacuum to an semi-infinite photonic crystal. 
As an example we recalculate a system studied previously by K. Sakoda et al. and get results in 
agreement with theirs. 

PACS numbers: 42.70.Qs,85.60Bt 



Since Yablonovitch 0, 0] discovered a periodic di- 
electric structure that has an absolute gap in the fre- 
quency spectrum for electromagnetic waves, the idea of 
photonic crystals have attracted great interest. Many 
phenomena have been predicted theoretically and many 
application possibilities have been explored EL IE 01- 
Corresponding studies in the early years most authors 
paid their attention on the frequency spectrum gaps, and 
the most convenient method to calculate the band gaps 
is the plane wave expansion method 0,0- Recently, var- 
ious kinds new methods have been proposed to compute 
some other relevant physical parameters such as trans- 
mittance and penetration depth [i^ Hcl ITll IT^ | for a finite 
system. 

In this paper we also address the transmittance and 
penetration depth problems, but use a different method 
that includes all the information of propagating and 
evanescent modes getting from the translation operator. 
We show that by appropriately modifying the conven- 
tional plane wave expansion method we can enlarge its 
application region, and which makes it easy to solve the 
interface problem between vacuum and a semi-infinite 
photonic crystal system. 

Our method has several advantages. First, both the 
"air rods in dielectric" and "dielectric rods in air" prob- 
lems can be solved, without any restriction on the shapes 
of the rods and position of cutting plane that separat- 
ing the semi-infinite photonic crystal region from the 
air region, which is impossible for the LKKR method 
HH [To . ITtI llgj . Second, all information getting from 
the complete set of the eigenfunctions of the translation 
operator are used, including both the propagating and 
evanescent modes. This makes it easy to analyze and 
discuss phenomena using the well established knowledge 
of solid state physics. Third, the finite size effects such as 



the resonance behavior of the transmittance curve caused 
by the finite thickness of the photonic crystal sample can 
be easily isolated. We can thus accurately calculate the 
transmittance and reflectance for a very thick photonic 
crystal sample. 

For a system without free charge and current, and if 
the permittivity e(r) and permeability /i(r) are scalars 
independent of time, the magnetic field H(r, t) satisfies 
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V x -V x H(r,t) 



(i) 



where H(r,i) = ]TH w (r)e- Mt , and 
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V x -V x HL(r) = ^H^r). 



(2) 



In addition, if e(r) and ^((r) are periodic functions, fol- 
lowing the derivation of Bloch theory, Eq. J2J can be 
changed to 

- (k + G) x e^_ G , (k + G') x Hk,c 

G' 

= u 2 ^fi G -G'H k ,G', (3) 

G' 

where H w (r) = Ee'' k+G > r H k , G , e (r) = £e iGr £ G , 
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fi(r) — E e' Gr /iG, and {G} is the set of the recipro- 

G 

cal lattice. Since in this paper we consider only two- 
dimensional cases, we have k z — 0, and the electromag- 
netic waves can be decoupled as E polarization (TE) and 
H polarization (TM) modes. For example, the TM mode 
of H field is written as H = H z z and satisfy 



'Electronic address: 
t Electronic address: 



ychsu. ep87g@nctu.edu. tw 



yangtj@cc. nctu.edu.tw| 



£>G-G' (k + G)-(k + G')if*,k,G< 
G' 

= ^ /iG-G'-ffz.k.G'- 
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Region I 




FIG. 1: A schematic view of the light incident from region 
I to region II, where region I is vacuum and region II is PC. 
The gray, dotted, and black arrows represent the incident, 
reflected, and transmitted light, respectively. For the incident 
light, k = {k x ,k y ). According to Bloch theory, the modes 
of region II can be written as H z {r) = J2g H G e %{ ^' +G) r . 
Based on the continuity conditions at the interface, reflection 
and transmission modes have k re / = ( — k" x ,k y + G y ) and 



ktrans = (k' x + G x ,k y + G v ), where k x = 
and k' x 's are obtained from solving Eq. 



— (ky + Gy) 2 



Conventionally, treating Eq. as an eigenvalue prob- 
lem, the propagating bulk modes of an infinite periodic 
system for a given real k can be obtained straightfor- 
wardly. However, in most of the situations we also need 
to know the transmittance and reflectance of a finite or 
semi- infinite system for a incident light. In order to ob- 
tain these quantities, various methods such like LKKR 
[H E E3, transfer matrix E E E3, and scat- 
tering matrix method have been proposed. In these 
kinds of methods, a photonic crystal slab is treated as the 
stack of many gratings. The matrix problem for only one 



grating layer is solved first, then multiplying the matrices 
layer by layer and the total transmittance and reflectance 
can be determined. Although these methods are success- 
ful, however, in order to confirm the good numerical ac- 
curacy, the number of layers should be restricted to a 
small value. In addition, using these methods it is hard 
to find the relations between the transmittance and the 
original band structure. 

In this paper, we propose an alternative method to 
calculate the transmitted and reflected waves from the 
interface between vacuum and a semi-infinite photonic 
crystal for a incident plane wave. The method is based 
on Eq. (0J , which contains all the information of the band 
structure of the system. We thus can easily analyze the 
physical meanings of various phenomena using the knowl- 
edge getting from the solid state physics. 

Here, if the k x and k y are fixed, the frequency u can 
be solved as an eigenvalue problem. If the system is in- 
finitely extended, there are just only propagation modes 
that can survive. However, sometimes we have to calcu- 
late the transmission and reflection coefficients for a fi- 
nite sized or semi-infinite sample, which have at least one 
boundary. At the edges of the sample, the periodic struc- 
ture is broken and the evanescent modes must be consid- 
ered. However, it is impossible to obtain the evanescent 
solutions from the eigenvalue problem of Eq. because 
it provides the solutions for an extended bulk and the 
boundary conditions at infinite restrict the k x and k y to 
be taken as real values. 

Nevertheless, can we make some modifications to pro- 
duce the attenuated solutions of Eq. (J2J? The answer is 
yes. For different purposes, there are two kinds of calcu- 
lations: one is to fix the direction of k and frequency; the 
other is to fix k y and frequency. With a simple transfor- 
mation, the Eq. Q can be rewritten as 
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for fixed direction and frequency, where (respectively) H G i, I and k denote the abbreviation of H z ^ G > ,5 G G > and 
unit vector of k, and 







£G-G" [w 2 /i G "-G' - e G '<-G' ( G " + k v9) ■ ( G ' + Ky)] P J V k x H G , 

I 



k x H G , 



(6) 



for fixed k v and frequency, where y and P denote unit 
vector of y direction and — £g-G" [ £ g''-G' 
respectively. If the k and G in Eq. J^J are acted by 
a rotation operator O which rotates k to x-direction — 
i.e. 6k = x — and we define G = 0G, then Eq. (0 



becomes Eq. (j^J in which G is substituted for G and 
k y = 0. Following this cue, the Eq. (0) can be consid- 
ered as a master equation for solving a problem where 
the incident light is always perpendicular to the inter- 
face. From which we can easily determine the penetra- 
tion depth along direction k. If our purpose is to use the 
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GaAs cylinder with r=0.1 5a,e =1 1 .43 




FIG. 2: The frequency spectrum of the square lattice pho- 
tonic crystal is derived from Eq. @. In this figure, solid lines 
and dash lines denote the TE and TM modes, respectively. 



band gap effect of the structure, then the result obtained 
from this calculation will tell us how to cut the sample 
to get the highest performance. 

On the other hand, Eq. © can be used to deal with the 
problem for light incident with different angle 0, which is 
the angle between the normal vector of the interface and 

k of the incident light. Here tan# = k y / .JuP/c 2 — k%. 

Figure n explains the details. 

Since we can obtain all the eigenvectors of the sys- 
tem, the transmission and reflection spectra can also be 
obtained. Based on the continuity conditions at the in- 
terface, the relationship between the H fields in region I 
and region II can be written as 



- (x y\Hl) (xov\H%) 
- (xovle-^Hl) (x y\e-^d x \H%) 



HL,\'R,\Hn 



(x y\Hl) 

(x y\e- 1 d x \H I ) / 



(7) 



where R and T are the reflection and transmission op- 
erators, and are the reflection modes in region 
I and transmission modes in region II, the m denotes 
the different modes, respectively. And Hq is the incident 
light. 

If Eq. J2J is expanded in K-space, it can be rewritten 

as 



i?i m |R|.Hio 
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where 



A = 
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-H 



I,r 

z,k.G„ 



(8) 



Ga-.G' 



II, t 
z,k,G 



and ff z ,k,G can be gotten from Eq. ©, and i,r and t 
denote incident, reflection and transmission, respectively. 
To determine the transmission and reflection coefficients, 
we have to first decide the direction of the Poynting vec- 
tor of every mode. For k is a real vector, it can either be 
obtained from v„ = VkW or from 



Re 



-HtVH z Ur 



cell 



G,G' 



Re <^ -H* 



z,k,G c G-G 



(k + G')#*,k,G' 



(9) 



For k is a complex number, propagation toward right 
hand side corresponds to Im(k x ) > 0. When group ve- 



locity and i? z ,k,G of each mode are known, the transmit- 
tance T and reflectance R can be gotten from them, and 
the accuracy can be estimated from how R + T is close 
to one. Sometimes, det |A| is possible to become zero, 
if so, the Eq. © has nonezero solutions when there is 
no incident light. This kind of wave is a surface state 
which resembles the surface plasmon propagating along 
the surface of a metal. However, we do not discuss it 
in this paper. One more interesting thing among these 
three equations — from Eq.(0} to Eq.© — is that they 
are identical to each other. Because the second row of 
matrices of left side of Eq. JSJ) and Eq. © are equal 
to Eq. (QJ while the {G} in these equations are equal. 
Thus any eigenfunction of one of these three equations 
satisfies another two equations. It leads to two useful 
things: (i) the real k contours of these three methods 
with equal frequency are the same, (ii) The {G} needn't 
to be changed when these three equations are considered 
as a series of policy tools. For example, we should de- 
cide where the band edge is and select the u near the 
band edge to obtain the band structure from Eq. if 
we hope to observe what happened near the band gap. 
By replacing the frequencies of Eq. (JSJ) and © by ujq, 
the penetration depth can be derived from Eq. JSJ and 
so do the variations of transmittance while the solutions 
of Eq. © are used in Eq. J7J. During this process, even 
if we just select G = and drop out G ^ for average 
assumption in the calculations of Eq. (0J, the {G} still 
needn't to be changed in the following calculations. 

For simplicity, the structure we used in obtaining Fig. [21 
and Fig. El is the square lattice with GaAs (e = 11.43 ) 
cylinders each with radius=0.15a inside vacuum; whereas 
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H-Polarization(real) H-Polarization(complex) 

(a) * - (b) 




im*(»))of E-Polarization 



FIG. 3: Possible k values for a constant-frequency lo — 0.4. 
The solid square frames in these figures are the first Brillouin 
zone boundaries and the data taken are inside the first Bril- 
louin zone, (a) and (b) are for the TM waves (H-polarization), 
(c) and (d) are for the TE waves (E-polarization). (a) and (c) 
are constant-frequency contours (with purely real k) of prop- 
agating waves, (b) and the inset of (d) show the contours 
of Re{k} of evanescent modes. In addition, (d) shows the 
min{Im £;(#)}. The spots in the (b) and the inset of (d) corre- 
spond to the cases with smallest Im(k) or longest penetration 
depth. 



in Fig. ^ we employ vacuum cylinders each with radius 
0.43077a inside the PbO(e = 2.72 ) background. In all 
cases, a is the lattice constant and the primitive vectors 
are given by ai = (a, 0) and a2 = (0, a). 

The first application is a selection rule to determine the 
interface direction for the highest performance of light 
insulation. We solve Eq. (J5J at a frequency equal to 
0.4(27rc/a) and the results are presented in Fig. [21 In 
Fig. |21 there is a band gap at lo — 0.4(27rc/a) for the TE 
mode, so that in Fig.^c), there are no real number so- 
lutions inside the first Brillouin zone. However, outside 
and far away from the 1st Brillouin zone such solutions 
exist, which are fake and are caused by the finite basis 
used in the calculations. To find the interface direction 
we first choose a direction k and use Eq. (0 to find a 
k that has the smallest |Im(fc)| value, which denoted as 
ki{9) and determines the main decay trend for a wave 
propagating along k. The second step is to scan angles 
from 0° to 45° to find an angle 6q that has the maxi- 
mum kj(6). The details are shown in Fig. |3fd), where 
we calculate the TE modes, and the penetration depths 
(i.e., 2n/ki{6)) for 0° and 45° are 6.8540a and 3.3272a, 
respectively. This indicates that when we produce a sam- 
ple that cut along 45°, it just needs 4 or 5 layers to stop 
the light with lo = 0.4(27rc/a) for TE modes instead of 7 
or 8 layers for 0°. 

The second application is to fix k y and frequency in 
order to obtain k x . For comparison, we select the system 
discussed in [j| to contrast with our system and show the 
result in Fig.0| The system in is a 16 layers photonic 
crystal which is a square lattice (with a lattice constant 



a) of air columns (radius equals 0.43077a and is located 
at the center of a unit cell) in a dielectric substrate PbO 
placed in air and our system is a semi-infinite photonic 
crystal placed in air. By using Eq. @ and Eq. © the 
transmittance can be obtained as shown in Fig. 01 Be- 
cause our system is infinite, we can find something quite 
different. 

(1) The solid lines and dashed lines are almost smooth 
curves except in the gap regime and at some special 
points (lo = 0.74 and 0.85 in TE mode). The oscillat- 
ing solid lines with dots in Fig. 4 represent the solutions 
of It is obvious that our curves are different from 
that of 0- The oscillation behavior is owing to finite 
thickness. They can be easily explained by a roughly 
consideration of the average dielectric e = (e) ce ;/- For 
low frequency, the most important contribution of £g is 
£g=0: which is equal to e, therefore the 16 layers photonic 
crystals can be considered as an effective material with 
uniform dielectric e = 1.7173eo and width 16a, where a 
is the lattice constant of the photonic crystal and eo is 
the permittivity of free space. This is a typical ID prob- 
lem and the waves in both sides of this material can be 
connected by transfer matrix whose formula is 



ikl6a 



o 





-ik!6a 



m: 



where M w is a 2 x 2 matrix with det |M W | ^ determined 
by the optical impedance contrast and the incident an- 
gle, and k = LO^/e/iQ is the effective wave number, where 
fiQ is permeability of free space. It is obvious that the 
transfer matrices equals ±1 when kl6a — nn, where n is 
an integer number. That means the frequency difference 
between two neighboring peaks is 



1 



1 



32 v /f7eo 32V1-7173 



0.0238, 



where lo ~ jp-. Comparing with the average width 
0.0222 of peaks between = 0.5 to 0.7, we can say 
the oscillation in Fig. 4 comes from the effect of finite 
size. 

(2) In the vicinity of Q ~ 0.74, according to the band 
structure calculation results the spectrum should ascend 
rapidly when lo is increasing, because it is at band edge. 
But, the real situation appears in Fig. 4(a) is ascending 
quickly and descending immediately to near zero trans- 
mittance. Our explanation is that there are two propa- 
gation mode's with group velocity v g ~ 0, so they do not 
contribute to the transmittance. 

(3) In Fig. 01a), the valley of transmittance near 
lo = 0.85 disappears when the interface is chosen to pass 
through the center of the vacuum cylinders. This re- 
veals that it is possible to stop the light at some isolated 
frequency points by appropriately choosing the cutting 
plane of the photonic crystal even if the frequencies are 
outside of band gaps. 

(4) In Fig.^Jb), there should be a forbidden band for 
lo = 0.75 to 0.78, but the line with dots does not show 
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0.5 0.6 0.7 

ra(27cc/a) 



FIG. 4: Transmittance of the (a)TE and (b)TM modes. The 
line with dots is the data excerpted from Ref. The solid 
lines and dashed lines are our solutions with different cut 
planes (in the solid line case the cut plane is the same as 
in whereas in the dashed line case the cut plane passes 
through the centers of the cylindrical holes), and the inset is 
the clearer view of (a) whose frequency regime is indicated by 



tion such like for how large a separation distance between 
two defects can they be treated as independent in super- 
cell method. 

In conclusion, the method we present here may not be 
efficient enough, because in the calculation we get results 
both inside and outside of the first Brillouin zone (FBZ). 
However, only the results inside the FBZ are useful and 
the others repeat the same information and are redun- 
dant. For example, if we use iV 2 bases, there are only 
27V useful eigenfunctions. Besides, this method has sev- 
eral advantages. First, from this method we can easily 
realize and analyze some properties of periodic systems 
with interface and the computational time is independent 
of the number of layers. Thus, even if the number of lay- 
ers is very large, it will save much time. Second, using 
Eq. @, we can also calculate the density of states, D(a>), 
through 



D{to) 



dkn 
|Vkw| 



(10) 



slie.ll 



this result. According our calculation the wave attenua- 
tion rate for a 16-layer structure is about 0.0733, which 
agrees with the result shown by the line with dots. This 
phenomenon shows that the evanescent modes do con- 
tribute the transmittance in a finite thickness structure. 

In these two applications, the evanescent modes are 
necessary and useful for the calculations of semi-infinite 
system, and this method also provides us some informa- 



They are especially useful when we aim to calculate the 
density of states in some small frequency regimes. 

We are now investigating the cases of finite size speci- 
mens and a structure with line defects. 

Finally, we thank Prof. B. Y. Gu for instructing us 
about Andreev reflection, which gave us a chance to em- 
ploy this idea about Eq. (JSJ, and thank Dr. P. G. Luan 
who let us find more possibilities with this method. 
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